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ABSTRACT 

We perform three-dimensional self-gravitating radiative transfer simulations of protoplanet 
migration in circumstellar discs to explore the impact upon migration of the radial temperature 
profiles in these discs. We model protoplanets with masses ranging between 10-100 M^, in 
discs with surface density profiles of 2 oc r"'^^, and temperature profiles of the form T cc r"^, 
where /3 ranges 0-2. We find that steep (fi > 1) temperature profiles lead to outward migration 
of low mass protoplanets in interstellar grain opacity discs, but in more optically thin discs the 
migration is always inwards. The trend in migration rates with changing /3 obtained from our 
models shows good agreement with those obtained using recent analytic descriptions which 
include consideration of the co-orbital torques and their saturation. We find that switching 
between two models of the protoplanet, one in which accretion acts by evacuating gas and 
one in which gas piles up on a surface to form an atmosphere, leads to a small shift in the 
migration rates. If comparing these models in discs with conditions which lead to a marginally 
inward migration, the small shift can lead to outward migration. However, the direction and 
speed of migration is dominated by disc conditions rather than by the specific prescription 
used to model the flow near the protoplanet. 

Key words: planets and satellites: formation - radiative transfer - methods; numerical - 
hydrodynamics - planetary systems: formation 



1 INTRODUCTION 

Investigations into the interactions of circumstellar discs and em- 
bedded protoplanets were first conducted by |Goldreich & Tremaine| 
( |1980j| . They found that the exchange of angular momentum be- 
tween the two should lead to migration of the protoplanet (latterly 
known as Type I migration), but their work did not suggest in which 
direction this migration would be. The analysis of |Ward| jl986^ 
for a non-self gravitating, two-dimensional disc, showed that if 
the disc had a negative temperature gradient (i.e. temperature re- 
duces with increasing heliocentric radius) then an embedded proto- 
planet should migrate inwards. When Hot Jupiters were discovered 
( [Mayor & Queloz|1995[ l and the difficulties with their in situ forma- 
tion became apparent, this inward migration seemed to be an ideal 
mechanism by which to reconcile existing formation scenarios with 
these planets' small final orbital radii. However, the timescales of 
growth and migration did not at first compare favourably. Migra- 
tion timescales were so short that protoplanets should plummet into 
their stars before they were able to grow to any considerable mass. 
Only by reaching masses sufficient to open disc gaps can proto- 
planets move from fast Type I migration to slower Type II migra- 
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tion, which proceeds at the disc's viscous evolution timescale. Once 
slowed, the planets only have to survive until the disc gas dissipates, 
thought to occur at a stellar age of less than 6 Myrs ( Haisch, Lada,| 
|& Lada|2001^ , at which point migration due to planet-gas interac- 
tions necessarily ceases. 



Since 'Ward (19 86)1 further ana lytical descriptions (i.e. [Ward] 
|1997[ [Tanaka, Takeuchi, & Ward] |2002[ l, and numerical mod- 
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have continued to find fast inward rates in the Type I regime (< 
100 Me). However, these works have generally considered locally- 
isothermal conditions. A large number of works, both analytical 
and numerical, have now been devoted to exploring the impact of 
more complex thermodynamics upon planet migration. The first of 



such works was published by |Morohoshi & Tanakaj ( [2003[ l, who 
performed local calculations focussing on a low-mass planet's in- 
teraction with its parent disc. They found that radiative cooling led 
to a non-axisymmetric mass distribution about the planet. This re- 
sulted in an additional torque beyond the commonly considered dif- 
ferential Lindblad torques, and altered the migration rate relative to 
a purely isothermal calculation. jMenou & Goodman) ( |2004) found 
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that a planet's migration rate was sensitive to the temperature and 
opacity structure of the disc through which it travelled, and under 
certain conditions the migration rate could be very much slower 
than that achieved in an isothermal model. This was closely fol- 
lowed by I Jang-Condell & Sasselov] ( |2005| l who also identified a 
slow down in Type I migration rates upon the introduction of radia- 
tive transfer. 

It was Paardekooper & Mellema ( 2006 1 who first found cases 
of outward migration due to a coorbital torque in their three- 
dimensional radiative transfer models. They found that this coor- 
bital torque could be prevented from saturating if the radiation dif- 
fusion timescale was shorter than the libration timescale of gas 
in the horseshoe region. This enabled temperature asymmetries to 
be maintained beyond a single libration period. They also showed 
that reducing the opacity (further shortening the radiation diffu- 
sion timescale) of their disc could lead to results in agreement 
with previous isothermal calculations. This was further described 
by |Baruteau & Masset ( 2008), who make clear that the radiation 
diffusion timescale must be greater than the period of a single 
horseshoe orbit to avoid returning to an isothermal like migration, 
and that outward migration can be related to the disc's radial en- 
tropy gradient. Many numerical models have since also found ev- 
idence of outward migration in the Type I regime (Paardekooper 
& Me lem a'2008nPaardekooper & Papaloizou'2008 ; Kley & Crida 
2008[|Kley , Bitsch, & Klahr 2009; Aylifll'e & Bate 2010, Yamada 
& Inaba|201 1( 1. Periods of outward migration can help to increase 
the overall migration timescale of a forming planet embedded in a 
disc, as is required by synthesis models to explain the population 
of exoplanets that has been observed jida & Lin|2008[ [Mordasim^ 
lAlibert, & Benz|2009| l. 

More recently there has been an extension of analytical de- 
scriptions of the non-linear coorbital torque, called the horseshoe 
drag, to describe it in both its unsaturated and saturated states 
l!Paardekooper et al. |2010[ [Paardekooper, Baruteau, & Kley|2011| 
[Masset & Casoli 2010|l. These works give expressions for the to- 
tal torque acting on a planet due to both Lindblad torques and the 
coorbital component. As such they can describe planet migration 
for a large range of scenarios, taking into account the disc's density 
and temperature profiles, as well as its thermal diffusivity. Outward 
migration is fastest in discs with steep radial temperature profiles, 
becoming more marginal at typical gradients such as T oc r"' . 

This paper is intended to explore protoplanet migration in a 
series of discs with radial temperature gradients from r" to r^^. We 
conduct three dimensional global disc models, using smoothed par- 
ticle hydrodynamics (SPH), that include self-gravity and which use 
a planetary surface to allow modelling of gas flow to well within 
the Hill sphere and the self-consistent formation of an atmosphere. 
We conduct a few isothermal models to investigate whether or not 
the temperature profile alone has an impact on the migration rate, 
but otherwise our calculations all include radiative transfer using 
a flux-limited diffusion approximation. In Section |2] we describe 
our computational method, in Section [5] we explain how we ob- 
tained our results which are then presented in Section|4] Section|5] 
discusses these results, whilst a summary and our conclusions are 
given in Section|6] 



2 COMPUTATIONAL METHOD 

The calculations described herein have been performed using a 
three-dimensional SPH code. This SPH code has its origins in a ver- 
sion first developed by |BenzH1990|[Benz et al.|1990[ l but it has un- 



dergone substantial modification in subsequent years. Energy and 
entropy are conserved to timestepping accuracy by use of the vari- 
able smoothing length formalism of Springel & Hemquist| ( |2002| > 
and Monaghan^(.2002) w ith our specific implementation being de- 
scribed in |Price & Bate| ( |2007^ . Gravitational forces are calculated 
and neighbouring particles are found using a binary tree. Radiative 
transfer is modelled in the two temperature (gas, T^, and radiation, 
r,) flux-limited diffusion approximation using the method devel- 
oped by Whitehouse, ^ate, & Monaghan|(2005j and Whitehouse] 
|& Bate] 12006^ 7 Integration of the SPH equations is achieved us- 
ing a second-order Runge-Kutta-Fehlberg integrator with particles 
having individual timesteps (Bate_1995p . The code has been paral- 
lelised by M. Bate using OpenMP. 



2.1 Equation of state and radiative transfer 

We present a few calculations performed using a locally-isothermal 
equation of state, as well as many more calculations which in- 
clude radiative transfer. In locally-isothermal models the temper- 
ature of the gas in the disc remains as a function of radius through- 
out the calculations. For the radiation hydrodynamical calculations 
we use the ideal gas equation of state, p = pTgRg/fi where Rg is 
the gas constant, p is the density, Tg is the gas temperature, and /j 
is the mean molecular mass. The equation of state takes into ac- 
count the translational, rotational, and vibrational degrees of free- 
dom of molecular hydrogen (assuming a 3:1 mix of ortho- and 
para-hydrogen; see |Boley et al.|[2007^ . It also includes the disso- 
ciation of molecular hydrogen, and the ionisations of hydrogen and 
helium. The hydrogen and helium mass fractions are X = 0.70 
and Y = 0.28, respectively, whilst the contribution of metals to the 
equation of state is neglected. More details on the implementation 
of the equation of state can be found in |Whitehouse & Bate|p006l ). 

The radiative transfer in these calculations is performed us- 
ing the flux-limited diffusion approximation, as implemented by 
[Whitehouse et al.| ( |2005l ) and [Whitehouse & Bate1j2006> , in which 
work and artificial viscosity (including both bulk and shear compo- 
nents) increase the thermal energy of the gas. Work done on the ra- 
diation field increases the radiative energy which can be transported 
via flux-limited diffusion. The energy transfer between the gas and 
radiation fields is dependent upon their relative temperatures, the 
gas density, and the opacity, k. Energy is lost from the system by 
the radiation field into the vacuum surrounding the disc; this is per- 
formed numerically through the use of a radiation boundary. This 
boundary is positioned at a height from the disc midplane at which 
the radiation path through the gas from this height outwards has an 
optical depth of Top ~ 1 ; i.e. the radiation can be expected to reach 
the vacuum. As a result, two layers of particles, one above and one 
below the disc, are made to comprise the radiation boundary. These 
particles are compelled to follow the initial temperature profile of 
the disc, causing them to function as energy sinks which impose 
a minimum temperature to which the disc can cool (see [Ayliffe &[ 
[Bate|20I^ for a fuller description). 

The use of steep radial temperature profiles in some of our 
models in this work leads to very high temperatures towards the 
inner boundary of these circumstellar discs. This material forms an 
unrealistically hot annulus at the inner boundary, from which radi- 
ation diffuses out through the disc, changing its structure out to the 
environs of the embedded protoplanet. To remedy this situation, the 
disc out to 2 au is compelled to maintain its initial temperature pro- 
file throughout the simulation, which quashes any diffusion in this 
region. This rule was applied in all the calculations for consistency. 
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including the shallower temperature profile calculations where it 
has no discernible effect. 

The opacities used here are those of Poll ack et al.| ( |1983| > and 
I Alexander] ( [ 1 975 j (the IVa King model), with the former providing 
the grain opacities and the latter the gas opacities at temperatures 
beyond the grain sublimation point. In some calculations the grain 
opacity is reduced by a factor of 100 to emulate possible modifica- 
tion of the population due to agglomeration processes (see [Ayliffel 
|& Bate|2009l for further details). 

2.2 Disc models 

We model a protoplanetary disc with radial bounds of 0.1 - 3 rp 
(0.52 - 15.6 au), where rp is the initial orbital radius of our embed- 
ded protoplanets, taking a value of 5.2 au. The disc is represented 
by 2 million SPH particles, a number found to deliver satisfac- 
tory resolution (see Ayliffe & Bate|20T0} . This leads to smoothing 
lengths at the disc midplane at rp of HjA (away from the planet), 
where H is the disc scaleheight and equals 0.05 at rp. Of course, 
considerably better resolution is obtained directly surrounding the 
protoplanet where the densities become much higher. The surface 
density of the disc goes as E oc r"''-, where the value at rp is £p = 
75 g cm"^ to allow comparison with previous works (e.g. Lubow, 



[Seibert, & Artymowicz|[T999l [Bate et~aL][2003] l. At the centre of 
the disc is a fixed potential representing a 1 Mq star. We imple- 
ment an inner boundary for the disc that prevents gas flow onto the 
star, which due to the lack of a self-limiting mechanism would oth- 
erwise artificially rapidly drain material from the disc. The outer 
edge of the disc is bordered by ghost particles which represent the 
disc beyond 15.6 au, preventing the disc from shear spreading into a 
vacuum (for more details see |Aylifi^e & Bate|2010[ l. The initial discs 
were evolved in the absence of a planet until any transience result- 
ing from settling had dissipated, which required just over 4 orbits 
of the disc's outer edge. Self-gravity is included in all the calcula- 
tions described in this paper. The impact of including self-gravity 
was discussed in |Ayliffe~& Bate (2010), where for these relatively 
low mass discs it was found to make only a marginal difference to 
migration rates, slightly slowing inward migration. However, self- 
gravity is essential for building a self-consistent atmosphere within 
the Hill radius. 

We perform calculations with a number of radial temperature 
profiles described by T" oc r"'' where j8 = 0, 0.5, 1, 1.5, and 2, and 
Jp ~ 73 K. Due to the use of steeper temperature profiles in this 
work, which results in higher temperatures at smaller heliocentric 
radii, the method used to set the disc temperatures is different to 
that used in |Ayliffe & Bate| ( |2010| ). In this previous work the initial 
temperature distribution was set assuming that the heat capacity at 
constant volume was independent of temperature, such that T oc u, 
where ii is the gas's specific internal energy. This is a good approx- 
imation so long as the rotational and vibrational modes of H2 are 
not excited (i.e. 7 = 5/3), which is true for T < lOOK. It was then 
the distribution of u that was actually prescribed to follow an r"' 
profile, such that yS « 1 ; this case is shown by the dashed lines in 
the left panels of Fig. [T] The resulting temperature profile is similar 
to the desired r"' profile as can be seen by comparing the solid line 
in the lower left panel, for which T is set explicitly, with the dashed 
line that results from setting 11. The higher temperatures in some of 
the disc models considered in this work fall beyond the range in 
which this assumption can readily be made, as shown by the right 
hand panels of Fig. [T] As such, for these models we set the initial 
disc temperatures and boundary temperatures to follow the desired 
profile (the solid lines in lower panels of Fig.[TJ and then the appro- 
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Figure 1. Internal energy (u, top panels) and temperature (T, bottom panels) 
profiles along the midplane of unperturbed discs. The dashed lines show the 
profiles that are given when it is u that is established with the given radial 
dependence. Conversely, the solid lines show the profiles that are produced 
when it is T which is compelled to follow the stated profile. For the shallow 
profile given in the left panels, the temperature profile that results from set- 
ting M oc r"' is similar to that obtained by forcing T oc r"' . However, where 
a steeper profile is used, as in the right panels, the temperatures reached 
are sufficient to bring about substantial changes in the specific heat capac- 
ity of the gas (i.e. due to dissociation of molecular hydrogen, excitation of 
rotational and vibrational modes of molecular hydrogen) such that the pro- 
portionality of energy and temperature breaks down; we use heat capacities 
calculated using Boley et al.| (20Q7l. As such in this work it is the tem- 
perature profile that is set directly (T oc r"''), for which the corresponding 
energy is then found. 



priate internal energy for the gas is found (solid lines in top panels 
of Fig. [TJ using tabulated heat capacities based upon [Boley et al.| 
(j2007| which are used in our radiative transfer implementation. We 
note that at the planet's orbital radius the disc's temperature is such 
that 7 « 5/3. Therefore we adopt this value for 7 when making 
comparisons with analytic models. 

The viscosity within the disc is introduced in a parameterised 
form, developed for SPH by |Monaghan & Gingol d (1983), which 
conserves linear and angular momentum, and was modified to deal 
with high Mach number shocks by |Monaghan| ( [l992) . It is a func- 
tion of two parameters, the o-spH-term establishes a viscosity to 
damp subsonic velocity oscillations that may be produced in the 
wake of a shock front, whilst the /?spH-term prevents particle in- 
terpenetration in supersonic shocks. This viscosity is only applied 
when the gas is under compression, and should ideally only act near 
shock fronts. We use the switch developed by |Morris & Monaghan| 
l |1997| l to try to reduce the action of artificial viscosity where the 
cause is not a shock. This attributes an o-spH-value to each particle, 
which is modified based on the local pressure gradient, such that it 
increases to a maximum of 1 at a shock, and falls away rapidly to 
0.1 as the particle moves away; this significantly reduces unwanted 
viscous dissipation in the disc. 
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3 CALCULATIONS 

3.1 Planet model and migration measurements 

Into the protoplanetary discs we embed a protoplanetary core, or 
partially formed protoplanet. In all the calculations discussed here 
the protoplanet is represented by a point mass which is free to 
move, and in all but a few cases (discussed below) the point mass 
is enwrapped by a surface of radius 0.03 times the Hill radius (Rn) 
( |AylifFe & Bate|2009||2010) . Gas is able to pile up on this surface 
to build an atmosphere, and so represent gas accretion in a reason- 
ably realistic manner. As in our previous work, a smooth start to 
the calculations is achieved by embedding a protoplanet of radius 
Rp = 0.0 Irp which then shrinks exponentially to the desired radius 
during the first orbit of the protoplanet. 

The point mass is free to move through the disc, enabling us to 
measure its rate of migration. This measurement is made by using 
a linear fit to the last 25 orbits of orbital radius evolution, in what 
is typically a 50 orbit total evolution. The first 25 orbits are given 
up to ensure that any transient disruption to the disc caused by the 
sudden introduction of a planet have subsided, and to allow gas 
in the horseshoe region to reach a quasi-equilibrium state. For a 
10 Me planet the libration time is » 50 orbits, for a 333 Mg planet 
it is » 8 orbits, whilst for the a 33 planet it is 29 orbits. As 
such, our lowest mass protoplanet mass models have completed a 
single libration period at the end of the calculations. In Ayliffe &| 
|Bate| pOTO) l we presented our results using migration timescales, 
T, which are calculated as the time a planet of fixed mass would 
take to migrate from its starting orbital radius (5.2 AU in all of the 
models discussed) in to the central star; t = r-p/r. This has no real 
meaning if the protoplanet is migrating outwards, so in this paper 
we present results using simply the migration rate, r, in units of rp 
per orbital period at rp (which for rp = 5.2 au is 11.86 years). 

We further refine the presentation of our results by quantify- 
ing the uncertainties that arise due to our method of measuring the 
migration rates. To achieve this, in addition to fitting over the en- 
tirety of the last 25 orbits, we also divide this interval into 3 smaller 
overlapping time domains to which individual linear fits are made. 
For a migration trend that is well represented by a linear fit, these 
3 migration rates should all be similar to both one another and to 
the rate measured over the entire 25 orbits. For less linear trends 
the measured rates give an indication of the possible variation that 
could be obtained by measuring at different points in the orbital 
evolution. We use the maximum and minimum migration rates ob- 
tained from these shorter time domain fits to provide error bars for 
each migration rate that is presented in this work. 

In [^iffe & BatelpOTO) we noted a difference between the 
migration rates measured when modelling a protoplanet with the 
surface described above and an accreting point mass. To explore 
this further we conduct a small number of calculations in which the 
protoplanet is represented as a point mass which accretes gas that 
comes within race, where r^cc = 0.03/?h in each case. The particles 
representing this accreted gas are removed from the calculation and 
their properties are used to modify those of the point mass, ensuring 
conservation of mass, linear momentum, and angular momentum 
(see |Bate, Bonnell, & Price|1995> . 



as the discs we have modelled in this work. They also include terms 
that describe the impact on the torques of the disc's thermal diffu- 
sivity, which in our radiative transfer calculations is effectively var- 
ied through changes in the opacity. In order that we might compare 
our numerical models with M asset & Casoli[ s analytic expression, 
we calculate the thermal diffusivity in our disc at rp in the absence 
of a protoplanet. In the flux-limited diffusion approximation of ra- 
diative transfer used in this work the flux in optically thick regions 
is given by 



3^P 



-V£, 



(1) 



where c is the speed of light, x is the Rosseland mean opacity, p 
is the gas density, and E = aT* in which a is the radiation density 
constant, a = 4cr^^/c. This allows us to restate equation[T]as 



F = ^^vr„ 

3^P 



(2) 



|Masset & Casoli| give an expression for the heat flux in terms of a 
form of thermal diffusivity, k, as 



F = -kVT. 



(3) 



Comparing equations|2]and[3]we obtain an expression for k, which 
is related to a true thermal diffusivity, k, by 



y - \ kfi 
y R„P 



(4) 



which is similar to equation 33 in [Masset & Cas oirf2010), where 
7 is the adiabatic index. The difference between this equation and 
that of lMasset & CasolTl lies in their use of E in the denominator of 
equation|4]as they are describing a two-dimensional disc, whilst for 
three-dimensions we use p. Making use of the approximation that 
p = 1,/H gives our final equation 



1 16o-sbT^^jH 



7 



3xpRI. 



(5) 



where H has a value of 0.05 at rp regardless of the chosen disc 
temperature profile. With this equation, taking 7 = 5/3 (which is 
appropriate at rp), and values from our unperturbed initial disc we 
are able to calculate the thermal diffusivities in models of differing 
opacity, x- For the 100% interstellar grain opacity disc this gives 
K = 6.57 X lO'^cm^s"', and the 1% opacity disc has a value of 
K = 6.57 X 10"cm-s"'. These diffusivities correspond with char- 
acteristic disc cooling timescales of tens of orbital periods and less 
than an orbital period, respectively, where this timescale is linearly 
dependent upon the chosen diffusivity. These timescales indicate 
that in the higher diffusivity (lower opacity) calculations the im- 
posed boundary will dominate the temperature structure of the disc, 
but this is not the case in the lower diffusivity (higher opacity) mod- 
els. 

It may also be of note that in these three-dimensional models, 
the dominant direction for radiative diffusion within the horseshoe 
region is vertically. For the high opacity calculations, vertical dif- 
fusion is found to be more than three times as efficient as radial 
diffusion. 



3.2 Opacity to diffusivity 

[Masset & Casoli| ( |2010| l developed an analytic expression to cal- 
culate the torques acting on protoplanets embedded in discs with 
temperature profiles described by the relation T = T^irlrpY^ , such 



4 RESULTS 

4.1 Is the surface treatment important? 

First we address a comment made in |Ayliffe & Bate||2010^ in which 
we suggested that the modelling of migrating protoplanets using a 
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surface treatment rather than an accreting point mass might signifi- 
cantly alter their migration rates, even changing the direction. This 
was inspired by a single result for a 10 Me protoplanet in a disc 
with a temperature profile exponent of /? « 1 , in which the migra- 
tion was seen to change from outward to inward when a crude evac- 
uating sink particle was replaced with a protoplanet with a surface 
in a high opacity disc (shown in the top panel of Fig.|2j. To explore 
this further we have conducted a series of calculations spanning 
10 - 100 Me in a disc with a steeper temperature profile, fj = 2. 
It was envisaged that such a steep profile would lead to outward 
migration for a low mass protoplanet modelled with either of our 
treatments by increasing the significance of the coorbital torques, 
making their effects more apparent. The results of these calcula- 
tions are shown in the lower panel of Fig.[2] It can be seen that using 
1% interstellar grain opacities (plus symbols), even with this steep 
temperature profile, the protoplanet migration is predominantly in- 
wards, with the 10 Me case being the only exception in exhibiting 
an almost stalled outwards migration. However, at 100% interstel- 
lar grain opacities those protoplanets with masses < 33 Me are all 
migrating outwards, regardless of the chosen protoplanet treatment. 
Comparing the migration of these low mass protoplanets with those 
seen in the top panel of Fig.[2]illustrates the increased importance 
of the coorbital torques in the disc with a steeper temperature gra- 
dient. 

In our new calculations (bottom panel) the 10 Me and 20 Me 
protoplanet migration rates obtained with an evacuating sink (un- 
encircled, black) are more rapidly outwards or less rapidly inwards 
than their counterparts modelled using sinks with surfaces (encir- 
cled, red); the reason for this is unclear, though it is in agreement 
with the sole case shown in the top panel for our previous calcula- 
tions. However, at masses > 33 Me this pattern reverses regardless 
of the disc opacity. |Bate et aT] ( |2003^ showed that it was for such 
masses that a protoplanet might begin to form a non-negligible 
gap in the disc. The evacuation of a gap reduces the mass in the 
coorbital region and as a result reduces the magnitude of coorbital 
torques. In |Ayliffe & Bate| ( |2010| we found that between the two 
treatments of the protoplanet used here, it was the evacuating sink 
particles that developed clearer gaps as a result of drawing in mate- 
rial artificially fast relative to protoplanets modelled with surfaces. 
This suggests that protoplanets of > 33 Me modelled with evac- 
uating sinks migrate faster inwards (slower outwards) than there 
equivalents modelled with surfaces because they reduce the coor- 
bital torques, which act outwards, to a greater extent. For a 100 Me 
protoplanet the coorbital region is sufficiently depleted with both 
treatments to lead to inward migration in a high opacity disc. 
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As to our findings in |Ayliffe & Bate| ( |2010^ that changing the 
protoplanet treatment can reverse the direction of migration for a 
low mass protoplanet, we must clarify that this appears to be a re- 
sult of the marginal nature of the migration of a 10 Me protoplanet 
in a high opacity disc with /? « 1. The change in the forces acting 
when comparing a protoplanet modelled using an evacuating sink 
and one modelled using a surface, happened to be sufficient in this 
case to change the direction of migration, though more typically it 
is likely to give a different rate with the same direction (as seen in 
Fig.[2|. The essential point is that as a protoplanet's migration rate 
is found to be at least somewhat dependent upon the method by 
which its accretion is modelled, it is best that the adopted treatment 
be as realistic as possible. We therefore remain convinced that us- 
ing a surface to allow gas-pile-up as a model of accretion is the best 
method for these purposes. 
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Figure 2. The upper panel shows migration rates, which were presented 
as migration timescales in'A yliflfe & Bate] j2010 '), for protoplanets mod- 
elled using surfaces (symbols encircled, red) embedded in discs with /3 » 1 . 
These discs were of either 100% (asterisks) or 1% (plus symbols) interstel- 
lar grain opacity. For a 333 Me protoplanet changing the opacity has such a 
small impact upon migration that the two rates appear indiscemibly overlaid 
above. The lower panel shows migration rates for protoplanets embedded in 
a circumstellar disc for which (3 = 2, some modelled with evacuating sinks 
(un-encircled symbols, black) and others using sinks with surfaces (once 
again, encircled symbols, red). All subsequent figures include only proto- 
planets modelled using sinks with surfaces. 



4.2 Impacts of the temperature profile 

Our primary goal in this work was to explore the impact of the tem- 
perature profile in a circumstellar disc upon the migration of pro- 
toplanets therein. In particular we wished to test the predictions of 
recent analytic models, which by their nature are unable to include 
details of the gas flow near the protoplanet; Fig. |3] contains the re- 
sults of our simulations to this end. Shown are the migration rates 
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Figure 3. Migration rates of four different mass protoplanets; across and down tlie panels 10, 20, 33, and 50 (as marked). Each protoplanet was modelled 
in a range of circumstellar discs with different temperature profiles, T oc r"^, where the temperature profile exponent fi is marked on the x-axis. Models in 
which interstellar grain opacities are used (asterisks) lead to slower or reversed migration when compared with the 1 % interstellar grain opacity models (plus 
symbols). pyiasset & Casoli ( 2010 1 give an expression for the total torque acting upon a protoplanet which includes the dependency upon the disc's temperature 
profile and its thermal dift'usivity, as well as the planet's mass. The resulting migration rates are plotted in each panel, one for a diffusivity akin to our 100% 
opacity calculation (solid line), and one for the 1% opacity case (dot-dash line). Migration rates obtained from |Tanaka et al!]j2002j are independent of the 
temperature profile, and so appear in each panel as a horizontal dashed line, the rate of which is determined by the disc surface density and the protoplanet 
mass. 



for protoplanets with masses 10, 22, 33, and 50 M^, embedded in 
discs with both high and low opacities, and a range of values of 
p. A horizontal dashed line in each panel marks the migration rate 
obtained using jTanaka et al.] ( |2002[ l's analytic form, which depends 
on the disc surface density and protoplanet mass, but is independent 
of the value of fi. It can be seen that the low opacity models ( 1 % 
interstellar grain opacity, marked with plus symbols) display some 
scatter, but generally centre around this analytic description. All of 



the low opacity calculations lead to inward migration, as would be 
expected for locally-isothermal calculations. 

Also marked in Fig.[3]are analytic migration rates based on re- 
cent work by Masset & Casoli (2010). In their description there is a 
dependence on bothyS and on the thermal diffusivity of the disc. As 
such we plot lines to correspond with both opacities used in calcu- 
lating our results, where the diffusivities used are those calculated 
in Section [3^ note that high opacities correspond to low thermal 
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diffusivities and vice versa. In the case of the high opacity calcu- 
lations (100% interstellar grain opacity, marked with asterisks) it 
can be seen that as the steepness of the temperature profile is in- 
creased, the inward migration rates slow towards zero before the 
direction changes to outwards and the rates begin to increase. In 
each panel of Fig. [5] excepting the 10 Me case for which the uncer- 
tainties are greatest, it is apparent that the rate of this change is in 
good agreement with the gradient given by the low thermal diffu- 
sivity (high opacity) result of Masset & Casolif s analytic descrip- 
tion (solid line). The analytic form also gives reasonable agreement 
(given the coarse spacing of j8 in our results) as to the temperature 
profile at which a given mass protoplanet in a high opacity disc 
will transition from inwards to outwards migration. The formula of 
[Masset & Casoli| is also dependent upon the disc's viscosity, which 
for plotting here we have used a |Shakura & SunyaevHI973 ) a-value 
of 0.004, which is equivalent to the magnitude of the SPH viscosity 
in the vicinity of the protoplanet in our models. Note also that this 
analytic description is valid only up to masses for which the half- 
width of the horseshoe region, x,, is less than the disc thickness 
at Tp. Adopting the adiabatic approximation of IBaruteau & Mas- 



|set| ^2008^ for Xj gives a maximum mass of 40 M^, meaning that 
the 50 Me cases shown in Fig. |3] are slightly beyond the range of 
validity. 

Using [Masset & Casoli[ s analytic form with values that cor- 
respond to our low opacity models indicates a reversal of the mi- 
gration trend with increasing resulting in migration accelerat- 
ing inwards for steeper profiles. However, the scatter in migration 
rates measured from our low opacity calculations is large, and pro- 
hibits us from concluding anything other that that the magnitudes 
show reasonable agreement. True locally-isothermal models were 
also conducted to establish whether they revealed a trend with the 
disc temperature profile. The results of these models are shown in 
Fig.|4] where once again there appears to be no definite trend in the 
change of migration rates with increasing /3, and the variation in the 
rates is small, varying by less than a factor of 2 from the minimum 
to maximum rate. As such we are unable to say whether or not the 
migration trend at low opacities given by [Masset & CasoUlpOTO) 
is real. 

Considering the immutable properties of our disc model, we 
find that only protoplanets < 100 M^ may migrate outwards for 
high values of p. This can be seen in Fig. |5] where migration rates 
for a 100 M(B protoplanet are plotted, and it is clear that no combina- 
tion of opacity and yS is able to cause outward migration. Included in 
the figure are the same analytic descriptions used in Fig.[3| but here 
the [Masset & Casoli[j201 ) line is well beyond the mass regime 
for which it is valid. For such a high mass protoplanet, the evacu- 
ation of the coorbital region should effectively end the significant 
influence of the coorbital torques. As such these effects are overes- 
timated by the analytic model, which is consistent with the fact that 
our high opacity results all lie well below the solid line. Noticeably 
the analytic model, even with the overestimation of the coorbital 
torques which should favour outward migration, does not predict 
such migration for this high mass protoplanet (for p < 2.5) , sug- 
gesting it is not possible in discs such as those modelled here. 

4.3 Impact of absolute temperature 

In this paper, we have studied the impact of the disc's tempera- 
ture profile on migration rates. In all these calculations, the ab- 
solute disc temperature at the planet's orbital radius is the same. 
As discussed above we find reasonable agreement with the ana- 
lytic description of [Masset & Casoli[ However, using [Masset &| 
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Figure 4. Migration rates for a 33 Me protoplanet in true locally-isothermal 
discs with varying radial temperature profiles. There is a small range of 
migration rates, varying by less than a factor of two from the minimum to 
the maximum rate. No significant trend is evident. 



[Casoli[ s analytic description it is possible to examine the migra- 
tion rate's dependence upon the absolute temperature at r^; this is 
presented in Fig.|6] The migration rates, calculated assuming a disc 
with a diffusivity of k = 6.57 x lO'^cm's"' (equivalent to our 100% 
opacity cases), are shown in the left hand panels. For each of the 
four protoplanet masses plotted (nip = 10, 20, 33, 50 M^) there 
is a temperature (scaleheight) at which the outward migration rate 
achieves a maximum, with this temperature scaling as rrP^^^ (within 
the mass range for which the model is valid, < 40 Me). The right 
hand panels of Fig. |6] show the same cases for a higher diffusivity 
oi K = 6.57 X lO'^cnrs"' (equivalent to our 1% opacity cases), for 
which conditions there is no outward migration, regardless of the 
chosen absolute temperature of the disc, and there is no turnover in 
the rate of migration, though these rates still vary significantly with 
temperature. 

These analytic results imply that as well as the temperature 
profile, the chosen initial disc temperature is an important factor in 
determining the migration behaviour of the various low mass plan- 
ets. Temperatures in different discs around various types of star may 
be expected to exhibit temperatures throughout the range shown 
here of 10-230K at 5.2 au, suggesting further modelling should be 
carried out for a range of disc conditions to explore the utility of an- 
alytic descriptions like ^Masset & Casoli[p010| l, which if valid, can 
be applied more generally than the results of numerical modelling. 



5 DISCUSSION 

First we address the low opacity calculations that we have per- 
formed. The resulting migration rates from these calculations ap- 
pear scattered, showing no evident trend with changes in the disc 
temperature profile. This can be seen in Fig. |7] which gathers all 
the low opacity results shown in Fig.[3]into a single panel, as well 
as including the analytic rates for an isothermal disc calculated 
using [Tanaka et al.[ l [2002^ (dashed lines) and rates based upon 
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Figure 5. Migration rates of a 100 Me protoplanet in a range of circumstel- 
lar discs with different temperature profiles, T oc r"'', where the temperature 
profile exponent /} is marked on the x-axis. Models in which interstellar 
grain opacities are used (asterisks) generally lead to very slightly slower 
inward migration than their equivalent 1% interstellar grain opacity (plus 
symbols) models. However, the trend with /J, which was quite pronounced 
for the high opacity models of lower mass protoplanets, is more subdued 
in this case. It also fails to follow the trend predicted using [Masset & Ca^ 
|soliH2010) for a high opacity disc (solid line), though at 100 Me we are 
applying the analytic model well beyond its range of validity; the analytic 
trend for a low opacity disc is marked with a dot-dashed line. The change in 
disc opacity leads to a smaller change in the protoplanet migration rate for 
this relatively high mass protoplanet when compared with the lower mass 
protoplanets shown in Fig.|3] Rates from jTanaka et al.]j2002^ are shown as 
a horizontal dashed line. 



the recent work of |D'Angelo & Lubow| ( [2010| l (solid lines) for a 
locally-isothermal disc. At these low opacities we expect the mi- 
gration results from our models to be similar to those predicted 
with isothermal analytic expressions as we found in |AylifFe & Bate| 
l |2010^ . Indeed, despite the scatter, these migration results do re- 
veal a trend with increasing mass that is matched by |Tanaka et al.| 
l |2002^ 's values. However, D'Angelo & Lubow have found that in 
the locally-isothermal (rather than globally isothermal) regime pro- 
toplanet migration rates depend upon the parent disc's temperature 
protile, with steeper negative temperature profiles (higher values of 
j8) leading to more rapid inward migration. This trend is not seen in 
our low opacity calculation results within which, as mentioned pre- 
viously, it is not possible to establish a trend that is distinguishable 
from the numerical uncertainties; this is also the case in the locally- 
isothermal calculations shown in Fig. [4] However, th e variation in 
migration rates obtained from |D'Angelo & Lubow| ( |2010| l across 
the range of /3 considered in this work is relatively small, less than 
a factor of 2 for each of the masses shown in Fig.|7] Our rather large 
uncertainties prevent us from definitively ruling out the model. 

Returning to [D'Angelo & Lubow| ( [20Tol >'s analytic form for 
an isothermal disc, it is interesting to note that the acceleration of 
inward migration with increasing yS is in qualitative agreement with 
the description of |Masset & Casoill ( |2010| l discs with high thermal 
difFusivity. As a brief aside we compare these two analytic models 
with which we draw comparisons. Fig. [8] illustrates the migration 
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Figure 6. Migration rates of protoplanets in discs of differing absolute 
temperature, calculated using the analytic description of [Masset & Ca-] 
|soli|p010) . Rates calculated for conditions equivalent to our 100% opacity 
models are shown on the left hand side, and 1 % opacity on the right. The 
temperature profiles are T oc r"' and r^- in the top and bottom rows re- 
spectively. The four lines in each panel from left to right correspond to 
protoplanets with masses 10, 20, 33, and 50 Me, respectively. The vertical 
dashed line marks the temperature at Tp used in our numerical models. 



rates calculated using |D' Angelo & Lubowf s isothermal formula for 
a 33 Me protoplanet (dashed line), with the rates obtained using 
[Masset & Casolif s formula for a whole range of thermal diffusivi- 
ties also shown. It can be seen that at the highest difFusivity (thick 
red line labelled C) the gradient of the change of migration rate 
with increasing /? for the two descriptions is extremely similar, and 
they are displaced in magnitude by a factor of a; 5/3, equivalent to 
7, as was found in [Paardekooper & Papaloizou| ( |2008] l. This similar 
behaviour is to be expected, with a high thermal difFusivity leading 
to a disc that resembles a locally-isothermal model, and responds 
to different temperature profiles in the same manner. As seen ear- 
lier when considering our low opacity calculations, the manner in 
which the migration rate oF a protoplanet changes with the discs 
thermal diffusivity is dependent upon the temperature profile of the 
disc. For discs with yS > 0.5, as the thermal difFusivity is increased 
From its lowest value (thick purple line labelled A) the migration 
rate oF a protoplanet can be expected to initially increase in the out- 
ward direction, or slow down iF the migration is inwards. Beyond a 
certain difFusivity, shown by the thick blue line labelled B in Fig. [8] 
the trend is reversed with increasing difFusivity slowing outward 
migration, and accelerating inward migration. This behaviour is a 
result oF the edge term oF the Horsehoe drag in |Masset & CasoIi} s 
model. In Future it would be desirable to explore a broader opacity 
space with greater refinement using our numerical models to deter- 
mine whether we can capture the effects oF the edge terms oF the 
Horseshoe drag. 

In our high opacity (low thermal diffusivity) models we find 
that increasing /? From leads first to slower inward migration and 
then accelerating outward migration For protoplanets oF < 50 Me. 
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Figure 7. Migration rates of four different mass protoplanets through low 
opacity discs of differing temperature profiles. The masses are 10 (plus 
symbols), 20 Me (asterisks), 33 Mg (diamonds), and 50 M® (triangles). 
Also included are analytic predictions. Those of |Tanaka et al.H 2002i with- 
out a/? dependence, for a three-dimensional isothermal disc, appear as hor- 
izontal dashed lines. The solid lines show the rates predicted using the /? 
dependent expression devised by [D'Angelo & Lubow| j2010) for a three- 
dimensional locally-isothermal disc. The masses that these lines correspond 
to are marked on the plot above the corresponding pairs. 



This is in agreement with the analytic work discussed above. Sim- 
ilar outward migration has been found in the numerical models 
of Paardekooper & Papaloizou (2008), Baruteau & Masset (20081, 
[Kley et al.^ ( ,2009j , and Paardekooper & Papaloizou (2009j , and has 
been associated with a density enhancement ahead of the planet 
and a deficit behind it in the horsehoe region. In each of these 
models the feature has been clearly identified in two-dimensional 
modelling, but such a feature is likely to be much less apparent in 
three-dimensional models where gas is not constrained to move in 
a single plane. The lefthand panels of Fig.|9]present plots of Er''- 
averaged over 25 orbits for a disc of /? = 2 containing a 33 pro- 
toplanet. The leftmost panel is a 100% opacity calculation in which 
the protoplanet exhibits outward migration, whilst the neighbouring 
panel is for a 1% opacity calculation that leads to inward migration. 
Relative to the low opacity case there is a more substantial deficit 
of material trailing the protoplanet (negative (p) in the horseshoe re- 
gion of the high opacity case. The rightmost panel of Fig.|9]shows 
a ratio of the two surface density perturbation plots, and makes it 
clear that there is a greater density contrast between the leading and 
trailing positions in the horseshoe region of the high opacity out- 
wardly migrating case. This density structure is reminiscent of that 
found in the two-dimensional models, though much less clearly de- 
fined. 



6 SUMMARY 

We have conducted a series of calculations to explore the impact of 
disc radial temperature profiles on the rates of protoplanet migra- 
tion. The surface density profile remains unchanged in each calcu- 
lation, thus each temperature profile gives rise to a different entropy 
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Figure 8. Migration rates based on the total torques calculated using 
analytic expressions for a 33 Me protoplanet in a disc with properties 
equivalent to our models. The dashed line indicates the expected locally- 
isothermal result for such a protoplanet in a three-dimensional disc, ob- 
tained using 'D'Angelo & Lu bow|j2010) . The solid lines are for rates ob- 
tained using Masset & Casoli ( 2010^ , where the colours progressing through 
the rainbow correspond to differing diffusivities, with the highest shown in 
red (labelled C). Three labelled thicker lines identify the trend with diffusiv- 
ity, with the corresponding values given in the key. The locally-isothermal 
and highest dift'usivity cases share the same gradient, with magnitudes dif- 
fering by K 5/3. 



profile through the disc. In agreement with analytic predictions, we 
find that steeper temperature profiles, /3 > 1, (and so entropy gra- 
dients) are much more disposed to bring about outwards migration 
for low mass (< 50 Me) protoplanets. There remains a dependence 
on the disc's opacity, which determines how the disc and the em- 
bedded protoplanet are able to cool, which can change the discs 
structure and the resulting torques. In interstellar grain opacity (low 
thermal diffusivity) discs, low mass protoplanets are found to mi- 
grate outwards, whilst in discs with opacities of 1% this level we 
find no cases of outward migration. The dependence of migration 
rates upon the disc's temperature profile in high opacity discs is in 
agreement with the predictions of recent analytic work by |Masset| 
|& Casoli| ( [20T0l l. 

For the low opacity calculations we find migration rates 
broadly inline with the rates predicted by [Tanaka et al.| ( |2002[). 
These conditions can be represented using the formulae of [Mas- 1 
|set & Casoli|p010l > with a suitably high thermal diffusivity. Such 
an approach suggests that a given protoplanet should migrate more 
rapidly inwards in a disc with a more steeply negative temperature 
profile. We found this to qualitatively agree with recent work by 
[D'Angelo & Lubowl l |2010) who find a similar trend for a locally- 
isothermal disc. However, our efforts to numerically test these pre- 
dictions were hampered by large numerical uncertainties in the re- 
sults of our low opacity calculations, leaving us unable to determine 
the veracity of the analytic predictions. 

We also expanded our investigation into the importance of the 
model used to represent the protoplanet. Using a surface that en- 
ables accreted gas to pile up and form an atmosphere results in 
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Figure 9. The two left panels show the surface density normalised by the density gradient of the initial unperturbed disc. Respectively the panels show the 
100% and 1% opacity models with jS = 2 and a protoplanet mass of 33 Mg. The right hand panel shows the ratio of the two left panels, revealing how the 
structures differ between the optically thick and thin models. This makes it clear that relative to the low opacity (inward migrating) model the high opacity 
(outward migrating) model has a greater contrast between leading (positive (ji) and trailing (negative 0) densities, thus explaining the dilfering directions of 
migration. 



somewhat different migration rates to those obtained when using a 
more crude accreting sink particle model. For the lowest-mass pro- 
toplanets (< 20 M(j,) the change in migration rate with protoplanet 
treatment was found to be small, though for borderline cases be- 
tween inward and outward migration it could be sufficient to change 
the direction. Those protoplanets of > 33 Mj, using the less realis- 
tic accreting sink model are found to migrate more rapidly inwards 
as they more effectively clear the coorbital region of material than 
equivalent surface models, more rapidly diminishing the positive 
coorbital torques. In general the effect of changing the protoplanet 
treatment is small compared to the differences brought about by 
altering the disc properties. 
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